load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

  varname = "NNUCCD"

  plotform = "eps"
  plotname = "profile_NNUCCD"

  ;;.......................
  ;; read data 
  ;;.......................

  fna = "pc_LP_dt1200.nc"
  fnb = "pc_LP_dt600.nc"
  fnc = "pc_LP_dt300.nc"
  fnd = "pc_LP_dt150.nc"
  fne = "pc_LP_dt75.nc"

  fla = addfile(fna,"r")
  flb = addfile(fnb,"r")
  flc = addfile(fnc,"r")
  fld = addfile(fnd,"r")
  fle = addfile(fne,"r")

  xa = fla->NNUCCD
  xb = flb->NNUCCD
  xc = flc->NNUCCD
  xd = fld->NNUCCD
  xe = fle->NNUCCD
  xt = fla->T
  xp = fla->PS
  hyam = fla->hyam
  hybm = fla->hybm

  ;;xa = where(xa.gt.1.e-20,xa,-999) 
  ;;xb = where(xb.gt.1.e-20,xb,-999) 

  xa@_FillValue = -999 
  xb@_FillValue = -999 
  xc@_FillValue = -999 
  xd@_FillValue = -999 
  xe@_FillValue = -999 

  nt   = dimsizes(xa(:,0,0,0))
  nlev = dimsizes(xb(0,:,0,0))

  v1 = xa(:,:,0,0)
  v2 = xb(:,:,0,0)
  v3 = xc(:,:,0,0)
  v4 = xd(:,:,0,0)
  v5 = xe(:,:,0,0)
  vt = xt(:,:,0,0)
  vp = xp(:,0,0)

  rho = vt 
  pre = vt 

  do ik = 0,nlev-1 
     pre(:,ik) = (/ doubletofloat(hyam(ik)) + doubletofloat(hybm(ik)) * vp(:) /) 
  end do 

  r = 287.042 ;; J/kg/K.

  rho = pre / (r * vt) 

  ;;.......................
  ;; #/kg/s -> #/L/s 
  ;;.......................

  ;;v1 = v1 * rho / 1.e3 
  ;;v2 = v2 * rho / 1.e3 
  ;;v3 = v3 * rho / 1.e3 
  ;;v4 = v4 * rho / 1.e3 
  ;;v5 = v5 * rho / 1.e3 

;;  v1 = where(v1.gt.1.e-9,v1,-999)
;;  v1@_FillValue = -999
;;
;;  v2 = where(v2.gt.1.e-9,v2,-999)
;;  v2@_FillValue = -999

  ;;.......................
  ;; define new vars 
  ;;.......................

  vv = new((/5,nlev/),"float") 

  printVarSummary(v1) 

  vv(0,:) = dim_avg(v1(lev|:,time|:)) 
  vv(1,:) = dim_avg(v2(lev|:,time|:)) 
  vv(2,:) = dim_avg(v3(lev|:,time|:)) 
  vv(3,:) = dim_avg(v4(lev|:,time|:)) 
  vv(4,:) = dim_avg(v5(lev|:,time|:)) 

  ;;vv = where(vv.gt.1.e-9,vv,-999)
  ;;vv@_FillValue = -999

  ;; inital radius 20um 

  ;;vv = vv * 4 

  ;;.......................
  ;; plotting
  ;;.......................

  wks   = gsn_open_wks (plotform,plotname)                   ; open workstation
  
  res                   = True                       ; plot mods desired
  res@tiMainString      = "" 
  res@trYReverse        = True                       ; reverse Y-axis
  res@xyDashPatterns    = (/ 0, 15 /) 
  res@xyLineColors      = (/ "red", "blue", "green", "orange", "purple" /) 
  res@xyLineThicknesses = (/ 4.0, 4.0, 4.0, 4.0, 4.0/) 
 
  res@pmLegendDisplayMode    = "Always"              ; turn on legend
  res@pmLegendSide           = "Top"                 ; Change location of 
  res@pmLegendParallelPosF   = .70                   ; move units right
  res@pmLegendOrthogonalPosF = -0.25                 ; more neg = down
  res@pmLegendWidthF         = 0.12                  ; Change width and
  res@pmLegendHeightF        = 0.10                  ; height of legend.
  res@lgLabelFontHeightF     = .016                   ; change font height
  res@lgPerimOn              = False                 ; no box around
  res@xyExplicitLegendLabels = (/" dt=1200 "," dt=600 "," dt=300 "," dt=150 "," dt=75 "/) 

  res@mpShapeMode  = "FreeAspect"
  res@vpWidthF      = 0.4
  res@vpHeightF     = 0.6

  ;;res@trXLog  = False ;;True
  res@trXMinF = -20. ;;1.e-6
  res@trXMaxF =  50. ;;6.e-3
  ;;res@tmXBLabelAngleF = 45	; tilt the XB labels 45 degrees

  res@tiXAxisFontHeightF    = 0.012 
  res@tiYAxisFontHeightF    = 0.012 
  res@tmXLLabelFontHeightF = 0.012 
  res@tmYLLabelFontHeightF = 0.012 

  res@tiXAxisString   = "NNUCCD (#/kg/s)" 
  res@tiYAxisString   = "Pressure (hPa)" 
 
  plot  = gsn_csm_xy (wks,vv,v1&lev,res) ; create plot

end



